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SUMMARY 

A fresh approach is taken to the embarrassingly difficult problem of 
adequately modeling simple pure advection. An explicit conservative control- 
volume formulation makes use of a universal limiter for transient interpolation 
modeling of the advective transport equations. This ULTIMATE conservative dif- 
ference scheme is applied to unsteady, one-dimensional scalar pure advection 
at constant velocity, using three critical test profiles: an isolated sine- 

squared wave, a discontinuous step, and a semi-ellipse. The goal, of course, 
is to devise a single robust scheme which achieves sharp monotonic resolution 
of the step without corrupting the other profiles. The semi-ellipse is partic- 
ularly challenging because of its combination of sudden and gradual changes in 
gradient. The ULTIMATE strategy can be applied to explicit conservative 
schemes of any order of accuracy. Second-order schemes are unsatisfactory, 
showing steepening and clipping typical of currently popular so-called "high 
resolution" shock-capturing or TVD schemes. The ULTIMATE third-order upwind 
scheme is highly satisfactory for most flows of practical importance. Higher 
order methods give predictably better step resolution, although even-order 
schemes generate a (monotonic) waviness in the difficult semi-ellipse simula- 
tion. But little is to be gained above ULTIMATE fifth-order upwinding which 
gives results close to the ultimate one might hope for. 

INTRODUCTION 

In a landmark series of papers in the 1970's, Bram van Leer worked 
"Towards the Ultimate Conservative Difference Scheme" for computational fluid 
dynamics (refs. 1 to 5). This work spawned a body of literature in the present 


‘Work funded under Space Act Agreement C99066G; presently at Dept, of 
Mechanical Engineering, The University of Akron, Akron, Ohio 44325. 



decade involving advective modeling methods which are sometimes classified as 
"shock-capturing" schemes or, more recently, as "TVD" schemes, referring to the 
oscillation-suppression strategy of total-variation diminution (ref. 6). The 
ultimate conservative difference scheme for CFD has proven surprisingly elu- 
sive. There have been notable "successes" such as impressive demonstrations 
of nonosci 1 latory high resolution of steps and shock waves, for example; but 
progress has been unexpectedly slow. It seems that correction of one defect 
always introduces another, equally severe. Unphysical oscillations inherent 
in classical second-order methods were eliminated by switching to first-order 
upwinding; but this merely replaced unacceptable oscillations with (what was 
ultimately realized to be) unacceptable global artificial diffusion. By devis- 
ing methods with locally varying artificial diffusion (small in smooth regions, 
larger in sharply varying regions), it Is possible to achieve somewhat better 
resolution than global first-order upwinding without introducing spurious 
numerical oscillations. 

Some forms of shock-capturing (or TVD) schemes achieve their impressive 
results for step resolution by the use of locally varying positive artificial 
diffusion or viscosity (first-order upwinding) to suppress oscillations, com- 
bined with local negative viscosity (such as first-order downwinding) to arti- 
ficially compress or steepen the front. Unfortunately, this inherent negative 
diffusion is responsible for artificial steepening of (what should be) gentle 
gradients, as well, as will be demonstrated. Because of the concomitant flat- 
tening of local extrema (due to the local positive artificial diffusion), this 
defect has become known as "clipping," although the problem is initiated by 
the artificial steepening introduced to give high resolution to simulated 
fronts. In some cases (for example, a step-function followed by a ramp), 
inherent oscillations, rather than being suppressed, are converted into a 
series of small monotonic steps, a phenomenon known as "stair-casing." 

As will be clearly demonstrated here, currently popular forms of TVD 
schemes achieve monotonic (although not always sharp) resolution of steps at 
the expense of gross (albeit nonosci 1 latory) distortion of simple smooth pro- 
files. Thus, A. R. Mitchell's characterization of advective modeling as compu- 
tational dynamics' ultimate embarrassment (ref. 7) is still appropriate, given 
the current state of the art. The ultimate goal of the current research pro- 
gram is to inject some self-confidence (as opposed to self-satisfaction) into 
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computational fluid dynamics by developing a truly robust (i.e., universally 
applicable) framework on which further refinements can be constructed. 

The present paper is the first in a proposed series presenting a fresh 
approach to simulating the advective term, which, being an odd-order (first-) 
derivative, is the most significant aspect of CFD (ref. 8). An extremely sim- 
ple form of "monoton i zing" universal limiter is described which can be applied 
to explicit conservative difference schemes, with no constraints on the order 
of accuracy and resolution. The universal limiter banishes unphysical over- 
shoots and nonmonotonic oscillations without corrupting the expected accuracy 
of the underlying method. In particular, when used with (artificial-viscosity- 
free) third or higher order base schemes, it does not induce artificial com- 
pression (steepening) or clipping, typical of so-called "high resolution" 
(actually, second-order) shock-capturing methods. 

Conservative explicit advection schemes of arbitrarily high order can be 
composed from "transient interpolation modeling" (TIM); i.e., for pure advec- 
tion of a scalar <t> at constant vector velocity v, the exact transient solu- 
tion over time At is 

<J>( x ,At) = <J>(x-vAt,0) (1 ) 

where the accuracy (in both space and time) of the approximate numerical method 
is determined solely by (multidimensional) spatial interpolation at the earlier 
time-level. The ULTIMATE conservative difference scheme then consists of using 
the universal limiter (UL) for transient interpolation modeling (TIM) of the 
advective transport equations (ATE). 

Since accurate modeling of the advective term is one of the more challeng- 
ing aspects of CFD (in addition to nonlinearities, multidimensionality, etc.), 
attention will be focused, in this first article, on the superficially simple 
but embarrassingly difficult problem of unsteady one-dimensional pure advection 
of scalar profiles at constant velocity. Three critical test profiles are con- 
sidered: an isolated sine-squared wave, representing smooth functions with a 

continuously turning gradient; a step discontinuity in value; and a semi- 
ellipse, combining discontinuous and continuous changes in gradient. At 
first, several (unlimited) explicit schemes are surveyed, together with some 
popular shock-capturing methods which are explained in terms of the Normalized 
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Variable Diagram (NVD). The NVD is also used as the basis for development of 
the universal limiter, which is then explained in terms of unnormalized vari- 
ables and applied to explicit polynomial TIM methods of second through eighth 
order. For most (relatively smooth) flows, the cost-effective third-order 
upwind ULTIMATE scheme gives excellent practical results. Local higher order 
resolution can be automatically built in with almost negligible additional 
cost. 


SCALAR ADVECTION 


Unsteady, one-dimensional, pure advection of a scalar <j> at constant 
velocity u is described by 


§4l_ _ u 

9t u 3x 


( 2 ) 


This equation can be integrated in time over a time-step At and in space over 
a control volume (CV) at station 1 from -Ax/2 to +Ax/2, assuming a uniform 
grid. This gives the exact conservative difference scheme 

♦"*' a ♦? - c («v - * 1 ) <3 > 

where the bars represent spatial averages, and the asterisks time averages, the 
superscripts designate time-levels, c is the Courant number uAt/Ax, and right 
and left time-averaged face values are indicated. The CV spatial averages can 
be written in terms of central node values plus a deviation term 

|. = 4». + DEV (4) 


Now assume that node values are related by an exact equation of the form 

C - ♦? • -c 0v - ♦») <5> 

where the face values now include the effects of the deviation terms. The 
deviation terms themselves will then satisfy 


DEV n+1 - DEV n = -c(DEV r - DEV^) (6) 

given suitable definitions of the terms on the right. Thus, from equations (3), 
(4), and (6), 
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( 7 ) 


C - ♦? - A(*r - DEV r) - (*l - DEV *)] 

which can be put in the form of equation (5) provided the effective average 
face value at any face f is defined as 

4> f = 4>f - DEV f (8) 

In this interpretation, equation (5) rewritten as 

<t>J +1 = - c(<|» r - (9) 

can still be considered exact for constant c, given appropriate definitions 
of the effective face values. In practical simulations, the face values are 
approximated and the form of the equation is extended to variable advecting 
velocity (Including sign-reversals), i.e., 

C - ♦" - ( c r<V ' c lO ( ' 0> 

Transient Interpolation Modeling 

One method of generating advective algorithms in the form of equation (9) 
is based on equation (1), which in one dimension can be written, for constant 
u, 

<|)J +1 = <|>(x,At) = <|><x-uAt,0) = 4> n (x-uAt) (11) 

where <t> n ( x > can be considered to be a function of the normalized local coor- 
dinate, i- = (x - x.)/Ax, and given by 

4> n (x) = «|»J + f(V (12) 

The homogeneous spatial interpolation function f can in turn be written as 

f(§) = (13) 

where g might be a polynomial with coefficients depending on differences (of 
various orders) of local node variables. For a conservative scheme, these 
coefficients can always be written as follows 

g(E) = <A i + 1/2 - a 1 _ 1/2 ) + (B ui/ 2 ' B i-l/2^ + • • • <14> 
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or 

9<V = h 1 + i/ 2 ( 5> - hj_ 1/2 <5> (15) 

where can be obtained from h. +1 ^ 2 by the indicated 

index by 1 on all involved node variables. Then, by defining 

V-5> - ", ♦,«<«> 

and 

V-e> - h i-i/2 <5> 

Equations (11) to (17) can be combined into a form identical to equation (9) 

4>J +1 = - c [<t> r ( c ) - V c)] (18) 

Note that the effective face values are in general functions of Courant number 
and that conservation is guaranteed by equations (16) and (17). (Extension to 
variable velocity parallels eq. (10), with face values being functions of their 
respective local Courant numbers.) 

For example, for first-order upwinding, transient interpolation modeling 
is based on a linear function 

f<?) = s(<J>" - for c > 0 (19) 

= for c < 0 (20) 

Thus, for positive or negative c, <j> r can be written in terms of quantities 
centered at the right CV face 

V c) = I [(♦"+! + ♦") “ SGN(c)(*J +1 - 4>")] = say (21) 

and, of course, 4>^ is obtained by reducing all indexes by 1. Note that the 
first term in brackets is the sum of node values straddling the CV face, 
whereas the second term involves the difference of these values. 

Similarly, the well-known Lax-Wendroff scheme (second-order central TIM) 
corresponds to centered quadratic spatial interpolation 


reduction in 

(16) 

(17) 
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( 22 ) 


f<5> - 5 


(*?+! - *i_i) (*U1 ~ 2 *j * *1-0 


and the corresponding right face value is 


V c> 







(23) 


valid for both positive and negative c values. Comparing equations (21) 
and (23), it is seen that these two schemes differ by the Courant-number- 
dependence of the coefficient of the first-difference term. Both give exact 
point-to-point transfer when c = 1: <|> r = <t>", <|> a = ; and 

similarly when c = -1 . 


The QUICKEST scheme (third-order upwinding TIM), is based on upwind 
biassed cubic interpolation (ref. 9) 


♦J 3> - - TTir K - »■<« $ < 

2 

where 19 r is the (right-face-centered) sum of second-differences centered 
at i and i+1 


19 


2 

r 





,n ,n , n , n 
= * 1+2 ‘ * 1+1 ' *1 + * 1-1 


(25) 


3 

and S r is the difference of the above second-differences, i.e., the third- 
difference (centered at the right face) 




cj) 1 ? 

m+2 


- 3<|> 


i + 1 


3<j>i - 


’ 1-1 


(26) 


Note that - SGN(c)6^/2 is the upwind biassed second-difference, propor- 

tional to the upwind curvature. Conversely, for fourth-order central TIM 

3 

(centered quartic Interpolation), less weight is placed on S r 


♦ 


(4) 

r 






) 


(27) 
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Second-order upwinding, based on upwind shifted quadratic interpolation, 
has a similar form to QUICKEST, but with a different coefficient of the 
upwinded curvature term: 

*< 2U) = 4 >J; 2> - ( 2^ 2 - SGN(c) S 3 ) (28 

and Fromm's method (ref. 10) is simply the average of equations (23) and (28), 
thereby cutting the curvature coefficient in half. 


Fifth-order upwinding and sixth-order central follow a similar pattern 


.(5) , (4) (1 - c 2 )(4 - c 2 ) 

♦r = ♦r + 2 T 5 I 


( 2 ^J - SGN(c) 


and 


.( 6 ) .(4) (1 - c 2 )(4 - c 2 ) 

^r = ^r + rrn 


K - § v 5 ) 


(29) 


(30) 


respectively, where 2 #* is the sum of fourth-differences 


2 K ■ *H3 - 3 *?.2 * 2 *?.1 * 2 +? - 3 *?-l * *?-2 


(31) 


and 6 J. is the face-centered fifth-difference 


^5 1 n p 1 n 1 rtfi 1 a . n r . n , n 

S r = *1.3 - S *1.2 * l 0 i.l - ,0 *1 * 5 *i-l - *1-2 


(32) 


Finally, seventh-order upwinding and eighth-order central are summarized 


by 


(7) , (6) (1 - c 2 ) (4 - c 2 )(9 - c 2 ) 


d> = <t» 
T r i t' r 


7! 


(l9 S r - SGN(c) sj) (33) 


and 


* (8) - V 


(6) (1 - c 2 ) (4 - c 2 )(9 - c 2 ) 


7! 


( 2 ®r - I S ?) <34> 


respectively, where 


* ♦?♦« * 5 *W3 * 9 *?.2 - 5 *U 1 - * 9 <-1 - 5 *?-2 * *?-3 


and 


(35) 


S r = ^ i +4 - <3 + 21 * U 2 - 35 * U 1 + 35 *" - 21 *?-1 + 7 *?-2 “ *?-3 


(36) 
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Ml 

Note that for Nth-order central methods, the coefficient of is propor- 

tional to c/(N/2), as compared with SGN(c) for the related (N-l )th-order 
upwind scheme. 

Central or upwinded polynomial TIM methods of the above forms can, of 
course, be continued up to arbitrarily high order via straight-forward recur- 
sion formulas, based on Binomial coefficients or Pascal's triangle. For the 
methods considered in this study, figure 1 shows the one-dimensional control- 
volume stencils involved as the order is increased. The stencils are shown 
for variable, possibly reversing, velocities, corresponding to equation (10), 
even though the test problems used in this study involve only constant veloc- 
ity, positive to the right. For this reason, first-order upwinding, for 
example, involves the same three-point stencil as second-order central, corre- 
sponding to equations (21) and (23), respectively. QUICKEST, second-order 
upwinding, Fromm's method, and fourth-order central each involve a five-point 
stencil, when velocity reversals are allowed for. Four grid points are 
involved for each face, as seen from equations (25) and (26). This is the 
same stencil used by all (second-order) shock-capturing TVD schemes. Clearly, 
as the order of accuracy is increased, a wider stencil is necessitated; an Nth- 
order central scheme and its related (N-l )th-order upwind scheme require an 
(N+l)-point stencil in one dimension. 

Test Problems 

The following three test problems are selected on the basis of simplicity 
and ease of reproducibility, and are intended to represent basic characteris- 
tics of behavior that might be encountered in practice. A given numerical 
scheme is considered successful if it is able to simulate all three test 
problems to within some desired level of performance; if a scheme fails one or 
more of the tests, it is deemed unsatisfactory no matter how accurately it 
simulates any one of the other tests. Performance is judged on the basis of 
two basic criteria: (1) total absolute error (ABSERROR) 

N 

g = | c i| <37) 

1 = 1 
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where e Is the local error at each node 


e * ^computed “ ^exact (88 

and (2) the WAVINESS or total variation of error 

N 

- «l| < 39 

1*1 

In addition. It Is usually desirable to monitor monotonicity; l.e., does an 
Initially monotonic profile remain so? In fact, strict monotonicity 
maintenance Is a basic (in fact the only) constraint determining the nature of 
the universal limiter to be developed here. 


The first test profile follows that used by Sweby (ref. 11), an Isolated 
sine-squared wave of width 20Ax 

4>(t*0) = sin 2 ^ 2 ^j for 0 < x < 20Ax 


(40) 


0 otherwise 


This function represents a relatively smooth profile with a continuously 
turning gradient and a single local maximum. In order to simulate practical 
situations, it is important to run the test problems over the same prescribed 
distance in all cases, irrespective of time-step (Courant number) or initial 
profile shape. In the tests described here, for example, the exact solutions 
advance by 45 mesh-widths. Two Courant numbers are used: c = 0.05 (900 time 
steps), representing "small" At; and a "moderate" value of c = 0.5 (90 time 
steps). The 25-time-step simulation used by Sweby was not long enough to see 
significant differences between the methods studied, or to allow their gross 
deficiencies to develop (which are typically worse at small Courant numbers, 
for a fixed distance) . 


The second test profile is a unit step change in 4> over one mesh width. 
Initially, <j> = 0 everywhere to the right of a specified jump point; all other 
points, including the upstream boundary, are set at 1.0. The unit step profile 
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is more fundamental than the isolated "square-wave," or box, used in some pre- 
vious studies (e.g., Sweby used an isolated rectangular box of width 20Ax in 
addition to the isolated sine-squared function of the same width). The box 
profile, at best, merely gives twice as much information as the unit step; but 
for highly oscillatory methods, oscillations excited by the step-up Interfere 
with those due to the step-down, and the resulting complex wave-pattern is not 
as enlightening as that of the simple step simulation. The unit step is also a 
basic test of monotonicity, a fundamental aspect of advective modeling. A 
"good" step simulation is one which monotonical ly resolves the step in a 
"small" number of mesh widths - the smaller the (monotonic) "numerical width," 
the better the method (for this test!). 

The third test profile follows one used by Steven Zalesak (ref. 12) that 

he attributes to B.E. McDonald. It consists of a semi-ellipse of width 

2i u Ax, initially centered at i r , 
w c 

♦t<t-0> ■ V' - <t - V 2/ ’w ,f I 1 - 'c| < 'w 

= 0 otherwi se (41 ) 

This is a rigorous test in that an Initial (leading) step change in gradient is 
followed by a region of continuously changing gradient and finally by a trail- 
ing step. Methods which are oscillation-free in the simple step simulation may 
generate significant waviness just behind the leading step or just ahead of the 
trailing step. The test used here differs slightly from that used by Zalesak 
in being 20Ax wide to conform to Sweby* s sine-squared profile, rather than 
30Ax ; this does not have any significant qualitative effect on results. 

Figure 2 shows the three initial profiles and the grid used for all tests. 
For reference, the leading-edge of each profile is positioned at the same loca- 
tion; the initial value is 0 for i > 23 in each case, and nonzero for smaller 
values. Grid-points 1 and 2 are used for boundary-condition treatment. Numer- 
ical boundary conditions are very simple: <t» a = 1 for the step profile and 

<|> a = 0 for the other two tests. Higher order methods (above fourth) require 
the designation of pseudonode values conforming to the above pattern. In all 
cases, computation begins with grid-point 1. To normalize the x-domain to 
order unity, Ax is taken as 1/100; but this is irrelevant to the computation, 
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since the only parameter is the Courant number c. As mentioned, all tests are 
run so that the exact solution would translate a fixed distance, taken to be 
45Ax . Thus, the total number of update steps, N^, is related to the Courant 
number by 



Results for Polynomial TIM Methods 

Figures 3 to 12 show the results of simulating pure advection of the three 
test profiles for first-order upwinding, the Lax-Wendroff scheme, second-order 
upwinding, Fromm's method, QUICKEST, fourth-order central, fifth-order upwind- 
ing, sixth-order central, seventh-order upwinding, and eighth-order central, 
respectively. Each case is run at both a representative small Courant number 
(c = 0.05, N t = 900) and a moderate Courant number (c = 0.5, * 90). First- 

order upwinding is the only polynomial method which gives monotonic step reso- 
lution; but the width of the resolution is poor and the other profiles show the 
well-known effects of this method's inherently large global artificial numeri- 
cal diffusion, which is worse at smaller Courant numbers, as seen in figure 3. 
The Lax-Wendroff method (fig. 4) generates trailing oscillations (phase-lag 
dispersion) typical of central-difference methods; in these cases, phase lag is 
worse at smaller Courant numbers (compare figs. 4, 8, 10, and 12). Second- 
order upwinding gives rise to leading oscillations, as seen in figure 5. It 
was this observation that led to the idea of averaging the Lax-Wendroff method 
with second-order upwinding in an effort to cancel phase-lead and phase-lag, 
at least in some "average" sense (ref. 10). Fromm's zero-average-phase-error 
method indeed shows marked improvement (fig. 6), with much less sensitivity to 
Courant number, as seen. Fromm's method actually cancels the leading disper- 
sion term in the truncation error only at c = 0.5; in this case it is identical 
to third-order upwinding (QUICKEST), as seen in the second half of figure 7. 
QUICKEST gives markedly better performance at other Courant numbers, but the 
two methods are qualitatively very similar. Note that the smooth (sine- 
squared) function is particularly well modeled. Fourth-order central is again 
highly oscillatory, with large ABSERR0R and WAVINESS for each profile, as seen 
in figure 8. The fact that fourth-order central methods are substantially 
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inferior to third-order upwinding is apparently not well known, and one con- 
tinues to hear of researchers who switch to central fourth-order schemes after 
experiencing "difficulties" with second-order. 

With fifth-order upwinding (fig. 9), one begins to see a trend which con- 
tinues to higher order: simulation of the smooth profile is excellent; as 

order increases, the step rise is steeper but the odd-order overshoots are 
larger and the even-order methods continue to be highly oscillatory, albeit 
with shorter wavelength; with the even-order central methods, significant wavi- 
ness develops in the semi-ellipse simulation just behind the initial jump in 
gradient. These trends are seen by scanning across figures 7 to 12. Finally, 
note that central (necessarily even-order) methods are much more sensitive to 
Courant number. This is because the highest order term in the face-value 
expression is proportional to c rather than SGN(c); compare equations (23), 
(27), (30), and (34) with (21), (24), (29), and (33), respectively. All the 
methods considered here give exact point-to-point transfer at c = 1 , as seen 
from the formulas for <j> r ; i.e., 4> r * <t>" , 4> a « ancl = For 

the higher order methods, point-to-point transfer also occurs (over more than 
one mesh width) for larger integer values of c; however, in the absence of 
modeled physical diffusion, these methods are not all stable over a continuous 
range, except for 0 < c < 1 . 

NORMALIZED VARIABLE DIAGRAM 
Normalized Variables 

Figure 13 shows a one-dimensional control -volume with attention focussed 
on one face (in this case, the left). In determining the effective face value, 
4y, the most influential nodes are the two straddling the face and the next 
upstream node, the latter depending of course on the flow direction at the face 
in question, i.e., the sign of u^. These three node values can be labeled <|> D 
(downstream), 4>y (upstream), and <j> c (central), as shown. Note the differ- 
ence in definition of these nodes, depending on SGN(u f .). In terms of original 
variables, there are clearly a very large number of cases to consider: combi- 

nations of positive or negative u^, positive or negative <J>, and positive or 
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negative values of gradient and curvature. Variations in sign, flow direction, 
and scale can be normalized out by defining the normalized variable (at any 
point) as 


$ 


4 > - <1>U 


(43) 


Now, in a case when <j>^ is a function of » <t>[j » and Courant number, 

the normalized face value is only a function of its adjacent normalized 
upstream node value and c 

* f - * c ) (44) 
since the other normalized node values are constant: 

$[] = 0 and $[J = 1 (45) 

Equation (44) Includes first-order methods, second-order central and upwind 
schemes, and third-order upwinding, in addition to second-order (and third- 
order) shock-capturing algorithms. Higher order methods involve more distant 
nodes but the normalized will still depend most strongly on $2 • 


Second-Order Time Averaging 

In general, for transient-interpolation-modeling methods based on 
equation (11), the time accuracy of the resulting control-volume algorithm 
(eq. (18)), is the same order as the degree of the spatial interpolation used 
in equation (12). Many advection schemes in common use (including TVD schemes) 
are based on second-order time averaging, which can be made explicit by 
straight-forward integration. Specifically, for spatially second-order- 
accurate methods, the deviation term in equation (4) is formally neglected; 
thus, from equation (8), the effective face values are just the estimated time 
averages. If these time averages are based on locally advected face-value 
linear behavior spatially, the method is second-order accurate in time, as well. 
For example if <jy is the instantaneous face value, the time-averaged face 
value over At i s 
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(46) 


♦f 


* 



At 


nAt 


d 0 


$^(t) dx = <(>jl - 


UAt 

2 



where <j>J! is the estimated face value at time-level n. In most cases, the 
spatial gradient is estimated as 


{3$\ n _ (*f ~ *c) 

V3x7^ Ax/2 

so that 

*f " *f " c (*f - «£) 

or, in terms of normalized variables, 

<jy = (1 - c)<j>£ + c^ 


(47) 


(48) 


(49) 


where, now, the linear Courant number weighting is explicit and depends 
only on <j>£ 


♦f = fn (^) (50) 

Figure 14 shows the normalized variable diagram (NVD) , plotting character- 
istics of the form of equation (50), for 


(1) First-order upwinding (1U): 


rn 


(2) The Lax-Wendroff method (20: 



(3) Second-order upwinding (2U): 



4 ? 


and 


(51) 


(52) 


(53) 
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(4) Fromm's method (2F): 


♦? - \* *£ <54) 

Note that the latter three spatially second-order methods all pass through the 
point (0.5, 0.75); in fact, any method whose characteristic passes 

through this point with a finite slope is (at least) second-order accurate in 
space, since can then be written 

= \ (l + ♦£) " CF ( 1 " 2 *c) (55> 

where the curvature factor CF is a constant or, in the case of nonlinear 
schemes, a function of <j>£. This is more apparent in terms of unnormalized 
variables 

*f - K*0 + *c) - CF (*o - 2 *c * *u) (56) 

or 

- CF Ax 2 + ... (57) 

thus, deviating from the second-order-accurate linear interpolation by second 
or higher order terms, provided CF is finite. Note that first-order 
upwinding, equation (51), and 

(5) First-order downwinding (ID): 

$£ = 1 (58) 

cannot be written in the form of equation (55) for finite CF. 

The characteristics shown in figure 14 are in the form of equation (50), 
i.e., the estimated face value at time-level n as a function of <j>£ . It is 
also important to portray the time-averaged face value in the same way, (no 
superscript), according to the linear weighting, equation (49). Figures 15 to 
17 show the NVDs corresponding to . c) for the Lax-Wendroff , second- 

order upwinding and Fromm methods, respectively. The characteristics are shown 
for five different values of Courant number: c ■* 0, c = 1/4, 1/2, 3/4, and 1. 

The zero-value is distinguished by a dashed line, since this can only be 
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approached in an explicit calculation. Note that all these methods revert to 
when c = 1 , as seen from equation (49), giving exact point-to-point 
transfer as discussed previously. 

Nonlinear "Shock-Capturing" Schemes 

For lack of better categorizing terminology, a number of currently popular 
algorithms have become known as "shock-capturing" schemes. When applied to the 
pure advection problems studied here, these second-order methods (in space and 
time) can be portrayed either in the <jy NVD or in the fy. NVD, correspond- 
ing to equations (50) or (49), respectively. However, in contrast to the 
linear characteristics of figure 14 (or of figs. 15 to 17), these methods are 
distinguished by nonlinear characteristics, although they remain linear in c 
according to equation (49). All schemes revert to first-order upwinding out- 
side the monotonic range, i.e., for ^ < 0 or <t>£ > 1 . They all pass through 
the origin (0,0) and the point (1,1) in either NVD. For the <}>£ NVD, they all 
pass through (0.5, 0.75), as required for second-order methods. 

Figure 18 shows the so-called minimum-modulus (Minmod) method (ref. 13) 
which is seen to follow second-order upwinding for 0 < 4>£ < 0.5 and Lax- 

Wendroff for 0.5 < <j>|l < 1. Part (a) of the figure shows <j>£ as a function of 

, whereas part (b) shows <jy as the c-weighted average between <j>J! and 
<|> c , equation (49), again for discrete values of c. A related method used by 
Chakravarthy and Osher (ref. 14) consists of second-order upwinding combined 
with first-order downwinding at time-level n. This is shown in figure 19. 

The MUSCL scheme of van Leer (ref. 5) is shown in figure 20; it consists of 

Fromm's method in the "smooth" (small curvature) region near 0.5, with 

piecewise linear deviations to pass through (0,0) and (1,1). In terms of <jy , 
the prescription for this method is 

(1) A sufficient (although not necessary) "monotonic" limiter: 

= 2 $£ for 0 < \ (59) 

(2) Fromm's method in "smooth" regions: 

rn 1 rn . 1 . rn 3 / rn\ 

4*^ = 4 + 4 £ 4 > q S. 4 (60) 
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and 


(3) First-order downwinding: 


= 1 for | < $£ < 1 (61) 

with first-order upwinding elsewhere, as usual. Another scheme due to van Leer 
(ref. 2) can be described by replacing the piecewise linear characteristic of 
MUSCL with a curved line (a parabola) for <j>£ in the monotonic range; this has 
the same slope as MUSCL for <j>£ = 0 + and 1_, as seen in figure 21. For con- 
venience, this scheme will be referred to as van Leer's "curved-line advection 
method" (CLAM); in the Sweby diagram (see appendix) the characteristic forms 
portion of a hyperbola, thus the scheme is sometimes referred to as van Leer's 
"harmonic" advection method. A scheme developed by Roe (ref. 6), nicknamed 
"Superbee," is shown in figure 22. Again piecewise linear, the 4>" character- 
istic (as <j>£ increases from 0) consists of a portion of slope 2, a portion of 
slope 1/2 (following Lax-Wendroff ) , a portion of slope 3/2 (second-order upwind- 
ing), and a portion of zero slope (first-order downwinding). This is considered 
to be one of the most "compressive" of the second-order schemes with respect 
to its narrow step resolution, as will be seen. This idea can be taken to 
extremes, however; another formally second-order scheme, which might aptly be 
called "Super-C" (for "compressive") is shown in figure 23(a); and an extremely 
"compressive" limited first-order downwinding scheme, "Hyper-C," is shown in 
figure 23(b). Note that both of these figures involve the time-averaged normal- 
ized face value, <jy. The Super-C scheme requires <jy to be the smaller of the 
upper universal limiter boundary (developed in the next section) and Lax- 
Wendroff, i.e., the smaller of 


<t> f = — and 



c) + j(l + c)$J 


( 62 ) 


for 0 < <j>£ < 1/2; with limited second-order upwinding, i.e., the smaller of 



and <)>^ = j(3 - c)<t>£ 


(63) 


for 1/2 < <j>£ < 1 , with first-order upwinding elsewhere. Hyper-C (limited 
first-order downwinding) simply requires <jy to be the smaller of 
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and 1 


(64) 


c 

in the monotonic regime, with elsewhere, as usual. Finally, note 

that the Courant-number-dependence of Super-C and Hyper-C is no longer linear. 

Results for the Nonlinear Schemes 

Figures 24 through 30 show the results of simulating the three test 
problems with each of the nonlinear schemes just described; again, the two 
representative Courant numbers (0.05 and 0.5) are used. As seen in figure 24, 
although Minmod resolves the step monotonical ly, the structure is not at all 
sharp and the other profiles are rather diffusive, especially at smaller 
Courant numbers. In the Chakravarthy-Osher simulations, shown in figure 25, 
the leading-edge steepening effects of the first-order downwinding are seen, 
with concomitant clipping due to first-order upwinding; again, profiles are 
more diffusive at small c values. The two van Leer schemes, MUSCL and CLAM 
(figs. 26 and 27), give similar results to each other, as perhaps expected from 
the qualitative similarity of their NVDs (figs. 20 and 21). The MUSCL scheme 
in particular is one of the more successful of the well-known second-order non- 
linear schemes, considering overall performance for all three test problems. 
Once again, both MUSCL and CLAM deteriorate at small Courant-number values, 
due primarily to the unnecessarily restrictive TVD limiter, equation (59). 

Superbee (fig. 28) gives significantly sharper results for the step simu- 
lation at all c values. The smooth-function (sine-squared) simulation is 
slightly better than that of MUSCL; however, there is a degree of artificial 
steepening inherent in this method, as seen in the semi-ellipse computation. 

In the presence of rapid changes in gradient (large curvature) - near the 
leading and trailing edges of the profile - the scheme has a tendency to con- 
vert all gentle slopes into sharp steps followed by plateaus. This defect is 
purposely taken to extremes with Super-C (fig. 29) and Hyper-C (fig. 30). 
Clearly, in terms of step performance, Super-C supersedes Superbee and Hyper-C 
supersedes Super-C. Super-C also does an excellent job of simulating the 
sine-squared profile; however, the semi-ellipse results are rather bizarre, 
showing stair-casing for small c values and gross artificial steepening at 
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c = 0.5. Although Hyper-C gives virtually exact step simulation at all c 
values, the other profiles are totally corrupted. These schemes were included 
to show the futility of designing a method on the basis of a single criterion 
(in this case, sharp monotonic resolution of steps); consequently the almost 
incredible step simulation of these artificial-compression methods should not 
be used as a standard for judging the overall performance of truly robust 
schemes . 

In order to show the effect of Courant number over a wider range, fig- 
ure 31 gives a log-log plot of ABSERROR for the sine-squared profile (lower 
curves) and the semi-ellipse (upper curves) versus ABSERROR for the step simu- 
lation for Minmod, MUSCL, Superbee, and Super-C, with Courant number as a 
parameter ranging from 0.01 to 0.978, with points shown at values of 0.1, 0.5, 
and 0.9 on each curve. At very small c values, Minmod produces large errors 
for all test profiles. The other schemes' semi-ellipse errors are comparable, 
with step errors decreasing in the order: MUSCL, Superbee, Super-C. The sine- 

squared error follows the same order. Superbee and, in particular, Super-C, 
show much less sensitivity to Courant number than the other schemes. Although 
all methods' trajectories in this plane approach the origin (exact point-to- 
point transfer) as c 1 , the tendency is rather slow, and even at c = 0.9, 
the sine-squared error in particular is unacceptably large, compared with what 
can be achieved with higher order methods. 

THE UNIVERSAL LIMITER 

The normal ized-variable plane, with <j>£ as abscissa and <j>f (no super- 
scrip) as ordinate, can be used to construct a very simple diagram representing 
constraints on the effective time-averaged normalized face value so as to guar- 
antee maintenance of monotonic profiles, thereby suppressing extraneous over- 
shoots or nonmonotonic oscillations, but allowing arbitrarily high resolution 
depending on the formal order of accuracy (in both space and time) of the base 
method. As seen in the previous section, nonlinear second-order schemes (lin- 
ear in c) can be represented by a single curve in the (<j>£ , plane for any 
fixed value of the Courant number. This is true for the nonlinear third-order 
scheme as well (to be described), but the Courant-number dependence is then 
also nonlinear, as will be seen. All the nonlinear second-order schemes 
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(except Super-C) satisfy Sweby's sufficient TVD criteria (ref. 10) (translated 
into present notation) 


= ijy = ij>£ for < 0 


and 


♦? > 1 


(65) 


and 


♦2 < <t>f < min^J, l] for 0 < $£ < 1 (66) 

In particular, all nonlinear characteristics pass through (0,0) and (1,1). 

These criteria are in terms of ; they are thus limited to temporally 
second-order schemes linear in c, according to equation (49). 


To allow higher order accuracy, it is extremely important to work directly 
with the effective time-averaged normalized face value <jy, rather than . 

By imposing simple monotonicity-maintenance criteria, much less restrictive 
constraints - the universal limiter - are placed on the allowable face values. 
The computational strategy is then, in principle, as follows: (1) formulate 

<l>f by some desired high-order method; (2) compute the actual normalized <j>£ 
value and the intended normalized value, (3) if falls within the 
allowable universal limiter range, simply proceed; (4) if <(y lies outside 
this range, reset (limit) its value to that of the nearest constraint boundary 
at the given <j>£ value; (5) reconstruct the new <jy from ; (6) repeat the 
procedure for all other faces; (7) update according to equation (10). In 
actual practice, it is less expensive computationally to carry out this strat- 
egy directly in terms of unnormalized variables, as will be described in 
detai 1 . 


Monotonicity-Maintenance Criteria 

Figure 32 shows normalized node values with <j>£ in the monotonic range, 
0 < < 1 • As suggested by the cross-hatching, monotonic behavior requires 
necessary conditions on <jy : 


♦r < 4> f < 1 (67) 

and on the corresponding face value of the adjacent upstream CV face, <j> u : 

0 < * u < $|! (68) 
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Consider equation (9), written for <j>£ in normalized variables 

*r = ? c ■ c (^f ■ o (69> 

In order to maintain monotonicity, the new <j>£ value must be constrained by 


n+1 

'U 




(70) 


For pure advection at constant velocity, the right-hand inequality is less 
restrictive than 4>^ 2 . but the left-hand inequality results in 


MVl (*c - C) 


(71) 


rn+1 

4>U 




and 

=0, 

~n+l 

♦u 

i .e. , 

nonpositive, the worst-case 

condi tion 

~n 

!c 

■ c 

for 

0 < <>[1 < 1 

(72) 

1 (67) 

y < 1 

for 

0 < <f>£ < 1 

(73) 


constitutes the universal limiter in the monotonic range of <j>£ . For <j>£ < 0 
or > 1 , numerical experimentation has shown that the simple condition 

4> f = 4>£ for $£ < 0 or $£ > 1 (74) 

gives the most satisfactory overall performance; this, of course, is equivalent 
to first-order upwinding as used by other nonlinear (second-order) TVD schemes. 
It does not erode the accuracy of the overall scheme, which is determined 
solely by behavior in the smooth region, near <j>£ ■* 0.5. 

The universal limiter, equation (74) together with inequalities (72) 
and (73), is shown in diagrammatic form in figure 33; the Courant-number- 
dependent boundary, = <j>£/c, is shown dashed to stress the fact that its 
slope changes with different values of c. Note that for c 0, this boundary 
approaches the vertical axis; while for c = 1, it degenerates into <j> f E <j>£ 
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(everywhere), corresponding to exact point-to-point transfer as usual. For 
reference to previous work, the corresponding criteria in terms of Sweby's 
variables, r and cp, are given in the appendix. 

For clarity, the precise steps in applying the universal limiter to 
transient interpolation modeling of the advective transport equations are 
given, as follows, using normalized variables. 

ULTIMATE strategy (normalized variables) : 

(1) Designate upstream (U), downstream (D), and central (C) nodes on the 
basis of SGN(u^) for each face in turn. 

(2) Compute DEL = - <f>J ; if | DEL | < 10" 5 (say), set <t> f = 4>£ and 

proceed to the next face. 

(3) Otherwise, compute <j>£ = - <|>^^/dEL; if this is less than 0 or 

greater than 1, again set ^ = <j>£ and proceed. 

(4) If not, compute <j>^ = - <j>^/ DEL, where 4>^ is based on a desired 

high-order TIM method. 

(5) If 4>f < <j>£ , reset (the lower limit) = <f>£ ; if <jy > <j>£/c, reset 
(the upper limit) <jy = <j>|?/c; if <jy > 1, reset (the absolute upper limit) 

= 1 • 

(6) Reconstruct <$y = <jy • DEL + <j>£j (this is the face value used in the 
update algorithm). 

(7) After finding all face values in this way, explicitly update according 
to equation (10). 

In order to avoid divisions and multiplications involved in constructing 
normalized variables and reconstructing the unnormalized variables, it is 
better to work directly with the unnormalized variables. 
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ULTIMATE strategy (unnormalized variables) : 


(1) Designate upstream (U), downstream (D), and central (C) nodes on the 
basis of SGN(u f ) for each face in turn, and compute DEL = 4 >q - 4>y and 
ADEL = | DEL | for each face. 


(2) Compute ACURV = 


i n « ,n I n 
<t>n “ t-vr + fii 


for each face; if ACURV > ADEL 


(nonmonotonic), set <j>^ = <f>£ and proceed to the next face. 

(3) Compute the reference face value <j> RE p = <j>[} + - 4>y)^/c for each 


face 


(4) Set up some desired high order face value . 

(5) If DEL > 0, limit <J> f by <j>£ below and the smaller of <t> RE p 

<t>p above. 

(6) If DEL < 0, limit <j>p by <t>£ above and the larger of <J> RE p 

below. 


and 

and <|>p 


(7) Update according to equation (10). 

In order to get some feeling for the ULTIMATE strategy, figures 34 to 37 
show normalized variable diagrams, <jy = f(<j>£ ,c) for the universal limiter 
applied to Lax-Wendroff, second-order upwinding, Fromm's method, and the third- 
order QUICKEST scheme, respectively. Note similarities (and differences) 
between figures 35 and 19(b), and between figures 36 and 20(b). Also note 
qualitative similarities between the limited versions of Fromm's method and 
QUICKEST (figs. 36 and 37) respectively; they are identical at c= 0.5 (as 
well as at c = 1, of course). 

Results for the ULTIMATE Schemes 


Figures 38 to 46 show the results for the three standard test problems 
obtained by applying the universal limiter to the Lax-Wendroff scheme, second- 
order upwinding, Fromm's method, QUICKEST, fourth-order central, fifth-order 
upwinding, sixth-order central, seventh-order upwinding, and eighth-order 
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central, respectively, each for the two Courant numbers c = 0.05 and 0.5. 

Note the inadequacies of the limited Lax-Wendroff scheme (fig. 38) and limited 
second-order upwinding (fig. 39) similar in some respects to the nonlinear 
shock-capturing schemes. Especially note the poor performance at the smaller 
c value; the reversed asymmetry of the profiles corresponds to the respective 
asymmetry of these two methods' NVDs (figs. 34 and 35). 

The limited Fromm method (fig. 40) is clearly an improvement over the 
simple second-order schemes, and has qualitatively similar performance to 
limited QUICKEST, as expected from the similarity of their NVDs (figs. 36 
and 37). The limited QUICKEST scheme (fig. 41) gives results which are 
probably entirely adequate for most practical situations. The sine-function 
error, in particular, is now within tolerable bounds. Although the limited 
fourth-order method (fig. 42) gives lower ABSERROR for each profile, there is 
a clearly discernible increase in WAVINESS in the semi-ellipse simulation, 
especially at small Courant number values; this is a typical shortcoming of 
even-order methods. Overall, the 1 imi ted-QUICKEST scheme seems to be the best 
of the schemes using the five-point stencil of figure 1; it is certainly far 
superior to any of the second-order shock-capturing schemes (which involve the 
same stencil). The artificial waviness of the limited fourth-order method 
(which also uses this stencil) detracts from an otherwise excellent 
performance. 

If better step resolution is required, the limited fifth-order upwinding 
scheme (fig. 43) gives highly satisfactory results, although of course this 
requires a seven-point stencil in general (allowing for velocity reversals). 
The limited sixth-order central method (which uses the same seven-point sten- 
cil), figure 44, gives slightly better step resolution, but worse performance 
for the semi-ellipse (in terms of both ABSERROR and WAVINESS) and a slightly 
higher ABSERROR in the sine-squared simulation at smaller Courant numbers (due 
to slight asymmetric clipping, typical of even-order methods). 

The higher order schemes follow a predictable pattern, with better step 
resolution, and essentially exact smooth-function simulation, but with annoy- 
ing waviness in the challenging semi-ellipse problem - now even noticeable in 
the upwind schemes, but still much worse with the even-order central methods, 


25 



especially at small c values. Clearly, the ULTIMATE strategy could be con- 
tinued to arbitrarily high order, either with polynomial TIM schemes of the 
type considered here or with alternate forms of interpolation such as splines, 
for example. The only stipulation is that the base method must be explicit in 
determining the intended - step (4) of the ULTIMATE strategy. 

Once again, to see the effect of Courant number over the complete stable 
range, figure 47 shows ABSERROR of the sine-squared and semi-ellipse profiles 
plotted on a log-log scale against ABSERROR of the step for the ULTIMATE Fromm, 
QUICKEST, fifth- and seventh-order upwind schemes, with Courant number as a 
parameter ranging from 0.01 to 0.978, with points shown at c * 0.1, 0.5, and 
0.9 on each curve. For clarity, the ULTIMATE even-order schemes have been 
omitted. The limited version of Fromm's method is perhaps of academic interest 
(being slightly better than MUSCL in the region near c = 0.5); but since the 
limited QUICKEST third-order upwind scheme requires the same stencil and 
essentially the same number and type of computations, it is clearly a more 
attractive method. Figure 47 should be compared with figure 31 giving the 
corresponding results for second-order shock-capturing methods. The obvious 
global characteristic for the higher-order ULTIMATE schemes is their much 
lower error for the smooth-function simulation, due to their lack of artificial 
steepening and concomitant clipping. As expected, step-simulation error 
decreases monotoni cal ly with the order of the underlying base method. The 
semi-ellipse and step errors of the higher order ULTIMATE schemes are strongly 
correlated, whereas for the second-order artificial compression methods of 
figure 31, the semi-ellipse error is roughly the same for each scheme, again 
reflecting the artificial steepening and clipping of these methods. 

Simplified ULTIMATE QUICKEST Strategy 

Referring to figure 37, it is seen that in the range 

0.2 < $£ < 0.8 (75) 

the ULTIMATE QUICKEST scheme is in fact identical to the unconstrained QUICKEST 
scheme. Thus, considerable simplification can be made in the algorithm without 
any approximation or effect on results. Inequality (75) can be rewritten as 

| - 0. 5 | < 0.3 (76) 
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or, multiplying by 2, 


|1 - 2$£| < 0.6 (77) 

Rewriting in terms of unnormalized variables results in 

| CURV | < 0.6 | DEL | (78) 

where CURV is the upwind-biased second-difference 

CURV = <i>J - 2<t>2 + 4>J (79) 

and DEL is the normalization difference 

DEL = <j>[J - <t>J (80) 

Thus, if inequality (78) is satisfied, the unconstrained QUICKEST scheme can 
be used directly, with no need for testing of universal-limiter constraints. 

In any practical flow, this criterion will be satisfied in the overwhelming 
bulk of the flow domain, being violated (if at all) only at a small fraction 
of grid points near where sharp changes in gradient occur. 

Of course, if 


$£ < 0 or > 1 (81) 

all ULTIMATE schemes (of any order) will use (in terms of unnormalized 
variables) 


<j>f = <J>£ (82) 

Inequalities (81) are equivalent to 

| $J! - 0. 5 1 > 0.5 (83) 

or, equivalently 


| CURV | > | DEL | 


(84) 



Thus, the simplified ULTIMATE QUICKEST strategy is as follows: 


(1) Designate "D," "C," and "U" nodes based on SGN(u^) in the usual way, 
for each face. 

(2) Compute | DEL | and |CURV|. 

(3) If inequality (78) is satisfied, use the unconstrained and unnormal- 
ized QUICKEST face value, 

♦f - K*D ♦ +c) - ¥ (*D - ♦?) - s 0 - c2 ) CURV <85> 

(4) Otherwise, if inequality (84) is satisfied, use equation (82). 

(5) Otherwise, compute the limited QUICKEST face value according to the 
(unnormalized) ULTIMATE strategy. 

(6) Proceed to the next CM face. 

(7) Update in the usual way. 

Note that step (5) occurs only in the small ranges 

0 < ^ < 0.2 or 0.8 < <j>£ < 1 (86) 


Cost-Effective High-Resolution Hybrid Scheme 

The ULTIMATE QUICKEST scheme is a simple, robust algorithm using the same 
stencil as second-order shock-capturing schemes but with much better global 
accuracy. If higher resolution of near-discontinuities is deemed necessary, it 
is clearly possible to use higher order (upwind) schemes globally. However, 
since the need for higher resolution occurs only in small localized regions, a 
cost-effective strategy is to use the efficient simplified ULTIMATE QUICKEST 
scheme as a base method, automatically switching to a higher resolution 
ULTIMATE scheme only where needed. 
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Numerical experimentation has shown that the need for a higher resolution 

scheme can be determined by monitoring the absolute value of the unnormalized 

2 

average "curvature," , defined in equation (25). If this is less than a 

specified constant, the limited QUICKEST scheme is used; if the critical value 

is exceeded locally, the algorithm branches to a higher order scheme (typically 

2 

fifth- or seventh-order upwinding), returning to ULTIMATE QUICKEST wherever 
drops below the threshold. By judicious choice of the threshold constant, the 
algorithm will use the ULTIMATE (probably unlimited) QUICKEST algorithm almost 
everywhere (i.e., in smooth regions) and switch to the higher order scheme at 
just the right grid points to give the desired high resolution. Figures 48 
and 49 show results for the usual test-problems for a hybrid-3/5 and hybrid- 
3/7 strategy, respectively, using a value of 0.015 for the threshold constant. 
Node values Involved in the higher order component on either left or right CV 
faces (or both) are shown by black dots; open circles designate that the 
simpler limited third-order scheme is to be used on both faces of the respec- 
tive control volume. 


SOME ASPECTS OF GENERALIZATION 

Clearly in this initial paper outlining the priciples of the universal 
limiter, attention has been narrowly focussed on the idealized academic (yet 
still very challenging) problem of one-dimensional pure advection at constant 
velocity on a uniform grid. This was done purposely, of course, to identify 
the basic problems associated with the advection term, the modeling of which 
is by far the most difficult numerical aspect and major pacing item in the 
development of computational fluid dynamics. [It makes absolutely no sense, 
logically or practically, to simulate a flow using highly sophisticated multi- 
equation turbulence models, for example, with an advection method which essen- 
tially replaces modeled physical viscosity with inherent (or explicit) 
artificial viscosity throughout the bulk of the flow domain - but this is still 
the "state-of-the-art" in a disturbingly large (and growing) number of commer- 
cial and research codes, especially in heat transfer and related industries.] 
But, obviously, in order for the universal limiter to be of practical value in 
a general-purpose code, several generalizations will need to be made. Space in 
a single (already lengthy) article does not permit detailed expositions of 
such generalizations or verification using classical test-problems. This will 
be taken up in future papers - both by the author and presumably by other 
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researchers who may wish to extend and apply the theory in various ways. How- 
ever, some generalizations are fairly straight-forward and will be sketched in 
this section, without showing specific results. Other generalizations of a 
more obvious nature (diffusion, nonuniform grids) or a more controversial one 
(systems of nonlinear equations) will be briefly addressed in the closing 
section. 


Variable Advecting Velocity 


The ULTIMATE strategy is easily extended to unsteady one-dimensional 
advection, where the advecting velocity is a function of space and time. For 
simplicity, assume that the right-face Courant number is positive 

c r > 0 (87) 


and that in a local region, 4> is increasing monotoni cal ly 


*1-1 < *1 < *1+1 < *i+2 


( 88 ) 


The update algorithm is based on equation (10), repeated here for convenience 


L n+1 


= - c r 4> r 


(89) 


where the face Courant numbers are considered to be known quantities at time- 
level n. As usual, <j> is first estimated using some desired high order 
method. This is limited by the adjacent node values 


< 4 > r < *J +1 


(90) 


Now require, conservatively, for monotonicity maintenance 


This can be rewritten, using equation (89), as 

c r*r * c a*d + *i " *i-l 

and using a worst-case estimate of 4 ^, this becomes 


(91) 


(92) 


c <J> < 

r Y r - 


n 


V1-1 


.H A 0 

+ 4 >. - 4 ) 


1-1 


(93) 
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for > 0. (If < 0, it may not be appropriate to require persistence of 
monotonicity.) Of course, for constant velocity, inequality (93) is equivalent 
to inequality (72). One further condition is necessary in the variable- 
velocity case, guaranteeing 


, n+1 


*1 1 * 1+1 


(94) 


which results in a condition on 4> a , for control-volume (i), given by 


± c r*r + *7+1 ' *1 ( 

or, in terms of 4> r , viewed as the left face of CV(i+l), using a worst-case 
estimate for the "far-right" face value, 

C <t> < C (b 1 ? . + <b? 0 - d>? . ( 

r^r - rrn + 1 m+2 m + 1 

assuming c rr > 0. In the constant-velocity case, this is superseded by the 
right-hand inequality of (90). The equivalent restrictions for monotonic 
decreasing regions should be clear. Local extrema are treated in the usual 
way. 


Nonl i near Advection 


Consider one-dimensional nonlinear advection 


where f(<t>) 
form as 


9$_. 9f( 4>) _ a 
8t 3x " u 


(97) 


is monotonic increasing. This can be rewritten in control -volume 


C' - ♦? * x ( f a - f r) (98 > 

where X = At/Ax, f^ = and f f = f(<t> r ^. Again, for definiteness, assume 

that 


*7-1 < *7 < *7 + i < *7+2 (99) 

First, estimate <j> r by some desired high order method. As usual, this will be 
limited by interpolative monotonicity 
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( 100 ) 


+1 * *r i 


which is equivalent to 

Xf(*J) < Xf(* r ) < Xf(*f +) ) (101) 

Then, to assure inequality (91), 

Xf r < Xf^) + ^ ^ (102) 

which is the generalization of equation (93). The condition analogous to 
inequality (96) is superseded by the right-hand member of equation (101). 
Clearly if 

f(<t>) = u<(> (103) 

for constant u, the limiting conditions revert back to inequality (72). 

Multidimensional Algorithm 

Because of the control -volume formulation of the one-dimensional 
algorithm, it is a straight-forward procedure to extend the ULTIMATE strategy 
to two and three dimensions. For two dimensions, the explicit update step 
analogous to equation (10) is 



. n+1 
*1 = 

= 

- ( c r*r - hh) - ( c b*b ' c t*t) 

(104) 

where "bottom" 

and "top" 

CV 

faces have been introduced, in addition to 

"left" 

and "right." 

Because of 

stri 

ct conservation, 






(105) 

and 








4> t < 1 . j > = 4> b < i .3 + 1) 

(106) 


(and similarly with the Courant numbers). Cell-averaged source-terms can be 
added to the right-hand side of equation (104), if appropriate. 

Focusing attention on one CV face (say the left), the first step is to 
compute some explicit high (third or higher) order mul tidimensional estimate 
for <$>£. This is then limited in a manner similar to figure 33, where is 
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based on three node values in a direction normal to the face: the two strad- 

dling the face together with the next upstream-biased node in the normal direc- 
tion, determined by SGN(u n ), where u n is the normal component of velocity at 
the face in question. Thus, in general terms, the limiting step is "locally 
one-dimensional" in the normal direction for each face, even though the high- 
order base method is multidimensional. 

For example, figure 50 shows the nodes involved for the left face (for u^ 
and v^ as shown) using the author's uniformly third-order polynomial interpo- 
lation algorithm (UTOPIA) as the base method, with the limiter-nodes shown as 
black dots and the remaining nodes as open circles. Allowing for all velocity 
directions on all four faces results in a 13-point stencil (the same as that 
used for the second-order discrete biharmonic operator). Higher order schemes 
can be designed in a similar way. Making use of additional nodes in the normal 
direction for each face appears to be more effective than involving other nodes 
in the transverse direction. Extension to three dimensions follows an entirely 
similar procedure. Details and results of rigorous test problems will be 
presented in a subsequent paper. 

DISCUSSION AND FORECAST 

The universal limiter, portrayed in figure 33, represents an extremely 
simple way to produce explicit "monotonic" advection schemes of arbitrarily 
high order. Potential numerical oscillations are suppressed without corrupt- 
ing the expected resolution of the underlying scheme. It has been seen that 
second-order methods (including well-known shock-capturing or TVD schemes) are 
significantly inferior to the third-order ULTIMATE QUICKEST scheme, which uses 
the same five-point stencil. Better resolution of sharp changes in gradient 
can be achieved by using a higher order base method, and concomitantly wider 
computational stencil. However, higher (even) order central methods are sus- 
ceptible to a degree of waviness under some circumstances, as seen with the 
semi-ellipse simulation. In terms of overall performance, there seems little 
to be gained beyond the ULTIMATE fifth-order upwind scheme, or the highly cost- 
effective hybrid scheme which uses ULTIMATE QUICKEST in smooth regions and 
automatically switches to a higher order ULTIMATE upwind scheme in isolated 
regions of high curvature. More precise simulation of narrow local extrema can 
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be built into the algorithm by using a weighting strategy which relies more 
heavily on an unlimited higher order scheme in such regions, provided a suit- 
able monitor can be devised to discriminate between physical and numerically 
produced extrema. This now appears to be possible, using simple pattern- 
recognition techniques. 

Some insight into the step-resolution of various schemes can be gained by 
correlating portions of a step simulation with regions of the (<j>£ , plane. 
Figure 51(a) shows schematical ly a simulated step profile, inverted for more 
direct reference to the NVD. In the leading-edge (LE) region, <t>£ is slightly 
less than 1; in the trai 1 i ng-edge (TE) region, it is slightly greater than 0. 
This is reflected in figure 51(b). As seen in figure 51(c), sharper resolution 
corresponds to larger ^ values relative to the corresponding <j>£ values, in 
both regions. This correlation can be quite clearly seen, for example, with 
the Chakravarthy-Osher scheme (figs. 19 and 25) or ULTIMATE second-order 
upwinding (figs. 35 and 39) where the NVD characteristic is relatively low in 
the TE region corresponding to a "blunt" trailing edge step simulation, and 
relatively high in the LE region, giving a sharp leading edge. The reverse is 
true for the ULTIMATE Lax-Wendroff scheme, resulting in a blunt leading edge 
and a sharp trailing edge. The Minmod scheme selects the lower characteristic 
in each region, leading to rather diffuse performance; whereas Superbee (or 
more effectively, Super-C) relies on the higher characteristic in each region 
of the NVD, and this results in the combination of sharp leading and trailing 
edges with concomitant narrow resolution of the step. 


The artificially high values of <j>^ for the latter schemes is equivalent 
to artificial compression, or (to the extent that the NVD characteristic lies 
above a well-behaved scheme such as ULTIMATE QUICKEST) to negative artificial 
diffusion. For example, the limited Lax-Wendroff scheme exhibits artificial 
negative diffusion for <j>£ < 0.5; whereas for limited second-order upwinding 
and the Chakravarthy-Osher scheme, this appears in the leading edge region, 
where the NVD characteristic is based on first-order downwinding. This 
interpretation correlates with the observed trailing oscillations of the 
unlimited Lax-Wendroff scheme (fig. 4) and the leading oscillations of 
unlimited second-order upwinding (fig. 5). 
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The great advantage of the ULTIMATE strategy is that it can be used with 
explicit schemes of arbitrarily high order, giving concomitant resolution of 
sudden changes in gradient without artificial steepening and clipping. The 
limited polynomial schemes reviewed here already give excellent results (above 
second order), but clearly the ULTIMATE philosophy can be used with other more 
sophisticated and more appropriate forms of spatial interpolation. One such 
method currently under investigation involves exponential ("tension") splines 
applied locally; this technique offers possibilities of very high resolution 
on a compact stenci 1 . 

Extension to multidimensional flow is straightforward and, as mentioned, 
this will form the subject of future papers. Although not addressed in this 
paper (since attention has been focused on the critical problem of pure 
advection) it should be clear that the introduction of modeled physical dif- 
fusion to an order consistent with the advection scheme presents absolutely no 
problems. The inclusion of physical diffusion in fact puts less demands on the 
advective algorithm; in most cases of practical interest, the ULTIMATE QUICKEST 
scheme gives results which are graphically indistinguishable from the exact 
solution (when known). In a similar way, extension to spatially varying grids 
presents no problems. The higher order accuracy of the base method can be 
formally maintained (refs. 9 and 16); or, provided the expansion ratio does not 
exceed about 125 percent, the simpler uniform-grid formulas for third- and 
higher order methods can be used directly without reducing the practical order 
of accuracy - although the formal order is, of course, reduced (ref. 17). 

One of the most important questions, of course, concerns how well the 
ULTIMATE strategy generalizes to handle systems of nonlinear equations such as 
the Euler (or Navi er-Stokes) equations of gasdynamics, and the "shallow-water" 
and related equations of oceanography and meteorology. This is by no means a 
trivial generalization and will no doubt form the subject of future papers; 
however, it is fair to predict that ULTIMATE simulations will be at least as 
successful as the better shock-capturing schemes, since the underlying philoso- 
phy is very similar. In particular, because of natural physical compression, 
shock-wave or hydraulic jump (or atmospheric front) resolution can be expected 
to be narrower than the scalar step resolution studied here. And because of 


35 



the ability to use arbitrarily high order resolution, unsteady gasdynamic simu- 
lations can be expected to give extremely reliable results, even at very high 
Mach numbers and nonideal conditions. Because of the conservative control- 
volume formulation, multidimensional gasdynamic, and geophysical algorithms are 
expected to be similarly robust. The general philosophy here is to model the 
control -volume advective modes precisely and let the physics (as reflected in 
the governing equations) take care of the "acoustic" or wave modes. This idea 
seems not to work well with currently popular second-order flux-limited schemes 
because of these schemes' addictive reliance on locally varying (positive or 
negative) artificial viscosity, requiring precise knowledge of characteristic 
speeds and directions. But the principle has great potential when used with 
higher order ULTIMATE multidimensional algorithms, which, by design, require 
only advective velocities at the CV cell faces, as prescribed by the govern- 
ing equations. 

Finally, the ULTIMATE strategy can be applied to steady-flow multi- 
dimensional simulations (even though the acronym might be somewhat out of 
place). Quite simply, the steady-state limiter consists of ULTIMATE with the 
Courant number set to zero. This, of course, results in the extremely simple 
condition in the monotonic regime 

4> c < <|>f < 1 for 0 < 4> c < 1 (107) 

with 4>f = <j>£, or similar characteristic passing through (0,0) and (1,1), in 
the nonmonotonic range. In other words, for locally monotonic behavior in the 
direction normal to a given CV face, after a high-order multidimensional 
estimate of 4>f is made, the face value actually used is constrained simply to 
lie between the adjacent upstream and downstream node values: <|>u and <j>d, 

respectively. The author's multidimensional SHARP algorithm (ref. 15), simple 
high accuracy resolution program, is a third-order nonlinear scheme conforming 
to these requirements; it is based on an exponential upwinding or linear extrap- 
olation refinement of the multidimensional QUICK scheme (ref. 16). But, of 
course, the principle can be extended to arbitrarily high order using the uni- 
versal limiter for tight resolution and accuracy, thus giving ULTRA-SHARP simu- 
lation of steady multidimensional flows containing near-discontinuities. In a 
practical algorithm, it is important to construct a single-valued upper boundary 
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for the allowable NVD region near 
by requiring 


4> ♦ 0 . 
C + 


This is most conveniently achieved 


<ty < const <j> c for 0 < 4> c 


(108) 


in addition to inequality (107), with suitably large values of the slope con- 
stant (~10 or 100). The piecewise linear NVD characteristic of the steady two- 
dimensional SMART algorithm of Gaskell and Lau (ref. 18) is also of this form. 
Steady-state methods based on second- and third-order schemes (including TVD 
limiters) can be solved with straight-forward ADI scalar penta-diagonal matrix 
algorithms using Gaskel 1-and-Lau 1 s modified (nonlinear) curvature-factor tech- 
nique, similar to the curvature factor, appearing in equation (55). 

Alternatively, arbitrarily high order limited schemes can be solved by popular 
ADI tridiagonal methods using the author's downwind weighting factor technique 
(ref. 19), thus rendering the ULTRA-SHARP very-high-order multidimensional 
nonosci 1 latory steady-flow schemes immediately compatible with many of the 
general-purpose commercial and research steady-state elliptic solvers currently 
in use, which are typically based on low-accuracy blended combinations of 
first- and second-order advection methods. 
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As described here (in the appendix), Sweby works with a much more restrictive 
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"TVD region," conforming to current practice in operational schemes, although 
the less restrictive conditions would have been known to him, since they had 
appeared in a GAMM proceedings paper by Roe and Baines (ref. 24) in 1982 and 
in other proceedings papers by Roe (refs. 6 and 25) in 1983. The author's 
attention was directed toward the latter papers by S.T. Zalesak in November 
1987 and by P.L. Roe in May 1988. Interestingly enough, Roe's 1986 review 
article (ref. 26) makes no mention of the less restrictive limiter boundaries, 
but rather describes instead only some of the more restrictive limiters (such 
as Minmod, "CLAM," and Superbee), even though the earlier proceedings papers 
are referenced. 

Zalesak's IMACS paper (ref. 12) describes work by his colleague, John 
Lyon, at the Naval Research Laboratory, who appears to be the only other 
researcher to have used arbitrarily high (up to eighth) order advection schemes 
constrained by an equivalent of the universal limiter. Finally, P.H. Gaskell 
and A.K.C. Lau, at the University of Leeds, have (again independently, follow- 
ing from their steady-flow work (ref. 18)) developed an equivalent limiting 
strategy for unsteady flow applied to a quasi-third-order advection scheme. 

They report highly accurate simulations of standard shock-tube test problems 
without resorting to flux vector splitting, Riemann solvers, or any of the 
other complexities associated with currently popular forms of first/second- 
order shock-capturing schemes. This confirms the author's belief that charac- 
teristic decomposition is not a necessity, but rather a symptom of poorly 
designed advection schemes, based on variable artificial viscosity, which 
require detailed explicit knowledge of all characteristic velocities and ampli- 
tudes. If the advection (macroflux) terms appearing in the conservation equa- 
tions for mass momentum and energy are carefully modeled with suitably limited 
higher order methods and microflux terms (pressure tensor, etc.) treated appro- 
priately, the governing equations themselves will automatically generate the 
correct additional wave modes. The far-reaching implications of this, espe- 
cially for multidimensional flows, should be clear. 
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APPENDIX 


In a well-known paper (ref. 10), Sweby introduced sufficient conditions 
for schemes to have total-variation-diminishing properties, and portrayed these 
conditions in what has become known as the "Sweby diagram," a plot of the flux- 
limiter factor, cp, against the local gradient ratio, r. For constant Courant 
number, the flux-limiter factor is defined as 


cp = 


Of - *c) 
0f w - *S) 


(A. 1 ) 


where <t> f is the limited face value, <|>£ is the adjacent upstream node value, 
LN 

and <jy is the face value according to second-order central differencing 
(Lax-Wendroff ) 



(A. 2) 


where 4>^ is the adjacent downstream value, as usual. All of the schemes 
considered by Sweby are linear in Courant number in the form of equation (49), 
rewritten for c > 0, as 

<j> f = (1 - |c|)4>J + |c|<J>£ (A. 3) 

The Lax-Wendroff face value can also be written in this form 

4>f W = (l - |c|)«t>£ EN + | c |4>2 (A. 4) 


where <Jy EN is the linear interpolation 1 + <j>£). Substituting (A. 3) and 

(A. 4) into (A.l) results in an equivalent expression for the flux-limiter factor 


cp = 


0 ? - ♦£) 
W EN - *9 


which can be rearranged as 


4 > 


n 

f 



+ cp 




(A. 5) 


(A. 6) 


interpreted as first-order upwinding modified by adding a term proportional to 
the difference between linear interpolation and first-order upwinding. 
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In terms of the normalized variables as defined in the present paper, it 
is not difficult to show that 


(♦f ~ *3 (*f ~ *c) 

K 1 - <3 " - ' c i ) (' - *3 


Similarly, the gradient ratio used by Sweby is 

_ (*c ~ ♦u) *c 

" (*D - *3 " 0 - *3 


(A. 7) 


(A. 8) 


The so-called TVD region proposed by Sweby consists of first-order upwind- 
ing (cp = 0) for negative r, followed by 


0 < tp < 2r for 0 < r < 1 


(A. 9) 


and then limited by first-order downwinding at time-level n, = 1 in equa- 
tion (A. 7), giving 

0 < cp < 2 for r > 1 (A. 10) 

Note that equation (A. 9) can be interpreted as 

♦y = 2 $£ for 0 < $£ < \ (A. 1 1 ) 

Figure A.l compares Sweby 1 s TVD region (shaded) interpreted in the (<j>£ , 
plane and in the (r, cp) plane. For reference, Roe's Superbee scheme is shown, 
passing through the "second-order" point, (0.5, 0.75) in figure A. 1(a), and 
(1,1) in figure A. Kb). By contrast, figure A. 2 gives the "extended Sweby dia- 
gram" corresponding to the ULTIMATE strategy of figure 33; recall that that 
figure involved <jy rather than <{>£ , as used in fig. A. 1(a). Note that the 
sloping boundary OB in the extended Sweby diagram has a slope of 2/ | c | as 
compared with the slope of 2 for 0C in the original Sweby diagram. Also 
note that the upper boundary of the extended diagram increases in height as 
| c | -» 1 ; in fact 

‘•’B-rrR 
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whereas, in the original diagram figure A. 1(b), the upper boundary is indepen- 
dent of c: 


<P, 


= 2 


As the Courant number varies, the point B in figure A. 2 "slides" along the 
line <p B = 2(1 + r g ), where r g = |c|/(1 - |c|>. For reference, the Super-C 
scheme is shown in figure A. 2. Hyper-C simply follows the upper boundary: 
AOBA. As is immediately obvious, Sweby's TVD region is grossly over- 
restrictive, resulting in predicted (normalized) face values which are too 
small, especially for values slightly larger than 0 or slightly smaller 
than 1. This is one reason for the poor performance of currently popular 
second-order TVD schemes. 
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(d) 

FIGURE 1. - CONTROL-VOLUME STENCILS FOR (a) FIRST-ORDER UPWINDING 
AND SECOND-ORDER CENTRAL; (b) SECOND- AND THIRD-ORDER UPWINDING, 
FROMM'S METHOD AND FOURTH-ORDER CENTRAL; (c) FIFTH-ORDER UPWIND- 
ING AND SIXTH-ORDER CENTRAL; (d) SEVENTH- AND EIGHTH-ORDER 
SCHEMES. 
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FIGURE 5. - RESULTS FOR 
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FIGURE 7. - CONCLUDED. 









































































FOR FIRST-ORDER UPWINDING (1U). LAX-WENDROFF (20, SECOND-ORDER UP- 
WINDING (2U), FROMM'S METHOD (2F), AND FIRST-ORDER DOWNWINDING (ID). 
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LAX-WENDROFF METHOD; Of AS A FUNCTI 
c — ' 0, c = 0.25, 0.5, 0.75, AND 1. 
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FIGURE 20. - NVDs FOR VAN LEER'S MUSCL. 



FIGURE 21. - NVDs FOR VAN LEER'S CLAM. 
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(a) c = 0.05. 

FIGURE 24. - RESULTS FOR THE MINMOD SCHEME 
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FIGURE 31. - ABSERROR FOR THE SEMI-ELLIPSE PROFILE (UPPER CURVES) AND SINE- 
SQUARED PROFILE (LOWER CURVES) PLOTTED ON A LOG-LOG SCALE AGAINST ABSERROR 
FOR THE STEP, WITH COURANT NUMBER AS A PARAMETER RANGING FROM 0.01 TO 
0.978, WITH VALUES SHOWN AT 0.1, 0.5, AND 0.9. CURVES SHOW RESULTS FOR: 
(1) MINMOD, (2) MUSCL, (3) SUPERBEE, AND (4) SUPER-C. ARROWS SHOW DI- 
RECTION OF INCREASING COURANT NUMBER. 
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FIGURE 39. - CONCLUDED. 
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ULTIMATE SIXTH-ORDER CENTRAL 
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FIGURE 44. - RESULTS FOR THE ULTIMATE SIXTH-ORDER CENTRAL SCHEME 
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FIGURE 44. - CONCLUDED. 


































ABSERROR (SINE-SQUARED) ABSERROR (SEMI-ELLIPSE) 



ABSERROR (STEP) 

FIGURE i\7. - ABSERROR FOR THE SEMI-ELLIPSE PROFILE (UPPER CURVES) AND SINE- 
SQUARED PROFILE (LOWER CURVES) PLOTTED ON A LOG-LOG SCALE AGAINST ABSERROR 
FOR THE STEP, WITH COURANT NUMBER AS A PARAMETER RANGING FROM 0.01 TO 
0,098, WITH VALUES SHOWN AT 0.1, 0.5, AND 0.9. CURVES SHOW: (1) ULTIMATE 

FROMM, (2) ULTIMATE QUICKEST, (3) ULTIMATE FIFTH-ORDER UPWINDING, AND 
W ULTIMATE SEVENTH-ORDER UPWINDING. ARROWS SHOW DIRECTION OF INCREASING 
COURANT NUMBER. 
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FIGURE 48. - HYBRID ULTIMATE THIRD-/FIFTH-ORDER UPWINDING RESULTS. SOLID 
CIRCLES SHOW NODES FOR WHICH EITHER THE LEFT OR RIGHT (OR BOTH) FACE 
VALUES FOR THE NEXT TIME STEP ARE TO BE COMPUTED USING THE HIGHER-ORDER 
SCHEME. 





























FIGURE 50. - TWO-DIMENSIONAL CONTROL-VOLUME SHOWING NODE STENCIL INVOLVED IN 
ESTIMATING THE LEFT-FACE VALUE USING ULTIMATE UTOPIA. SOLID CIRCLES SHOW 
NODES USED FOR THE UNIVERSAL LIMITER. 



(a) DIFFUSE SOLUTION. 




(C) SHARP SIMULATION. 


FIGURE 51. - SCHEMATIC (INVERTED) STEP SIMULATION SHOWING LEADING-EDGE (LE) AND TRAILING- 
EDGE (TE) BEHAVIOR. 
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(b) THE (M0) PLANE. 



FIGURE A. 1. - SWEBY'S TVD REGION (SHADED) SHOWN. FOR REFERENCE# ROE'S SUPERBEE 
SCHEME IS SHOWN BY THE HEAVY PIECEWICE LINEAR CHARACTERISTIC. 



FIGURE A. 2. - "EXTENDED SWEBY DIAGRAM* CORRESPONDING TO THE ULTIMATE 
STRATEGY OF FIGURE 33. FOR REFERENCE, THE SUPER-C SCHEME IS SHOWN 
BY THE HEAVY PIECEWISE LINEAR CHARACTERISTIC. 
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